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We discuss wave propagation along a crystalline piezoelectric fiber 
composed of lithium niobate or some other material in the trigonal 3m 
crystal class. The crystalline c axis is aligned with the fiber axis. We 
obtain an analytical description of all the vibrational modes. The 
method used is to make perturbation expansions about the modes of 
a hexagonal 6mm piezoelectric fiber, for which exact solutions are 
known. 

I. INTRODUCTION 

A single crystal of lithium niobate, grown in the form of a long fiber, 
has been considered for use as a low-loss acoustic delay line. Lithium 
niobate is of special interest because it is piezoelectric: it becomes elec- 
trically polarized when strained and, conversely, becomes strained when 
placed in an electric field. This piezoelectricity provides a means for 
electrically generating and detecting acoustic signals. 

In this paper we study mathematically the vibrational properties of 
a LiNb(>3 crystal fiber, with the crystalline c axis aligned along the fiber 
axis. The problem is by no means simple. We illustrate this by giving a 
brief history of related problems for which exact solutions have been 
obtained. The elastic, or acoustic, wave equations for an infinitely long 
circularly cylindrical isotropic rod were solved exactly by Pochhammer 1 
in 1876 and independently by Chree 2 in 1889. Even for an isotropic 
medium, exact solutions for a rod of finite length have not been obtained. 
It was not until 1965 that the next full exact solution was found. This 
was done by Mirsky, 3 ' 4 who determined the vibrational modes of a cir- 
cularly cylindrical rod consisting of a nonpiezoelectric medium which 
is transversely isotropic. Such a medium belongs to the hexagonal system 
of crystals; the crystalline c axis was aligned along the rod or fiber. Cer- 
tain of the modes obtained by Mirsky, i.e., those which are azimuthally 
symmetric about the fiber axis, had also been obtained earlier. 5 - 6 Re- 
cently, the author and J. A. Morrison were able to solve the coupled 
acoustic and electromagnetic wave equations, in the customary quasi- 
static approximation, for piezoelectric transversely isotropic crystals 
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belonging to the hexagonal 6mm, 622, and 6 crystal classes. 7 Exact so- 
lutions were obtained for all the vibrational modes. The author is una- 
ware of any other exact solutions, either for other crystals, or for other 
orientations of transversely isotropic crystals. 

The difficulty lies in the acoustic wave equations which, for a general 
anisotropic medium, consist of three coupled wave equations for the 
three vector components of displacement. If piezoelectricity is added 
via the quasistatic approximation, for which the electric field is repre- 
sented by the gradient of a potential, there are four coupled equations 
for four unknown functions. The boundary conditions may also involve 
all four functions coupled together. For the general anisotropic case, no 
method has been discovered to decouple the equations. For the specific 
crystals and orientations discussed above, it was possible to express the 
elastic displacements (and electric potential) in terms of three (or four) 
potential functions for which the wave equations decoupled. 

Unfortunately, such a serendipitous situation does not exist for the 
lithium niobate fiber. It belongs to the trigonal 3m crystal class; we 
cannot expect to find an exact description of the vibrational modes. It 
will be possible, though, to find an approximate description by means 
of an infinite series perturbation expansion. We use a technique which 
' is an extension of one used by the author to describe waves travelling 
along a sapphire fiber. 8 Sapphire is a nonpiezoelectric material belonging 
to the trigonal 3m crystal class. It is characterized by a stiffness matrix 
(used in the stress-strain relations) which has almost the same form as 
that for a transversely isotropic material. There is one additional stiffness 
coefficient. Since it turns out to be small in magnitude compared to the 
other stiffness coefficients, it is possible to describe the vibrational modes 
of a sapphire fiber (with the crystalline c axis aligned with the fiber axis) 
by means of perturbation expansions about the modes of a transversely 
isotropic fiber. 

The situation for LiNbC>3 is similar, albeit somewhat more compli- 
cated. We will make an infinite series perturbation expansion about the 
known solutions for a hexagonal 6mm crystal. The same techniques, 
incidentally, can be used to describe vibrations of crystals in the trigonal 
32 classes. We restrict ourselves to a discussion of trigonal 3m crystals 
only to keep the analysis from appearing extraordinarily complicated. 

For the sapphire fiber, numerical results are available for the low- 
est-order torsional mode of vibration; they are presented in a paper by 
the author and M. A. Gatto. 9 A low-frequency asymptotic analysis for 
that mode was also performed by R. N. Thurston and the author. 10 Ex- 
cellent numerical agreement between the results of the two independent 
theories provides a check on the rather complicated analyses involved 
and encourages us to extend the perturbation technique to a study of 
LiNb0 3 . 
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In Section II we write down the basic equations of motion and 
boundary conditions. In Section III we apply the perturbation technique 
and introduce potential functions. In Section IV we solve the differential 
equations, and in Section V we sketch how to apply the boundary con- 
ditions. 

Although it would be desirable to present numerical results as well, 
we shall not do so. Numerical results are not yet available for the un- 
perturbed (hexagonal 6mm) problem. The computational effort required 
to describe quantitatively the vibrations of a lithium niobate fiber would 
be even greater than the considerable effort expended to present results 
for a sapphire fiber. 

II. FORMULATION 

Consider a single crystal of LiNb0 3 (or some other member of the 
trigonal 3m crystal class), grown in the form of a fiber of circular cross- 
section, with the crystallographic c axis along the fiber axis. We shall 
assume that the fiber is infinitely long and has radius R. We adopt a 
cylindrical coordinate system whose z axis coincides with the fiber 
axis. 

In the quasistatic approximation, where the rotational part of the 
electric field is neglected, the basic differential equations are 11 

V • D = 0, (2) 

where T is the stress, D is the electric displacement, u is the elastic dis- 
placement, and p is the density. The properties of the specific crystal are 
introduced by means of the constitutive relations 

T = -e • E + c:S, (3) 

D = € • E + e:S, (4) 

where 

E = -V*, (5) 

S = V s u. (6) 

Here E denotes the electric field, S the strain, and $ the electric poten- 
tial. The crystal is described by means of the elastic stiffness matrix c, 
the piezoelectric stress matrix e, and the dielectric permittivity at con- 
stant strain matrix e. For a crystal in the trigonal 3m class, these matrices 
have the following forms in cylindrical coordinates: 12 
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Let n denote a vector normal to the fiber surface, i.e., in the radial 
direction. For the three mechanical boundary conditions, 13 we shall 
specify either that the surface tractions vanish: 



T*n = at r = R (free surface), 
or that there is no displacement at the surface: 

u = at r = R (clamped surface). 



(12) 



(13) 



The free surface condition is the natural one to consider for an acoustic 
delay line; it is equally simple to show how to solve the problem for the 
clamped surface condition, so we include it, too. 
For the electrical boundary condition, 13 we take either 



or 



* = at r = R (short-circuit), 



D*n = at r = R (open-circuit). 



(14) 



(15) 



The problem is to solve the four differential equations (1) and (2), in 
conjunction with eqs. (3) to (11), subject to four boundary conditions 
chosen from (12) to (15). Since we are concerned with waves travelling 
down the fiber, we assume the solution has an exp [i(ut — fiz)] depen- 
dence, where a> is the angular frequency and /? is the propagation con- 
stant; P will depend upon co. 
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We begin by writing the differential equations and boundary condi- 
tions in dimensionless form. Let 

Cij = Cij/c, c = max \cij\, 

U 

eij = eij/e, e = max|e i; |, 

hi = f.yA, t = max (t xx , t zz ). (16) 

Normalize u with respect to R, <t> with respct to Re/e, /3 with respect to 
/? _1 , and uj with respect to (c/p) 1/2 /R. To simplify notation, we use the 
same symbols as we used for dimensional quantities, except for the hats 
on Cij, eij, and ey. Upon substituting eqs. (3) to (11) into (1) and (2), we 
can write the dimensionless differential equations in cylindrical coor- 
dinates as 

en (u rr + - u r - — u J + c 66 — uee + (w 2 - 2 c 44 )u 

- 2i0c u cos 3d-u - 2i'0Ci4 sin 30 lu r - - uj 

+ (C12 + c 66 ) - v r0 - (cn + c 66 ) — v 
r r* 

- 2if3c u cos 30 (v r v) + 2i(3c u sin 30 - u e 

- if3(c l3 + c 44 )w r + 2c 14 cos 30 (- w r0 - — w 9 J 

+ C14 Sin 38 (w rr -~Wr ^ U>eo) 

- i0T(e x5 + e zl )$ r ~ 2re y 2 cos 3d (- $ rfl - — <f> g J 

- Te y2 sin 36 ($ rr $ r - — $oo) = 0. 

(C12 + c 66 ) " u r0 + (c„ + c 66 ) — u„ - 2ipc u cos 30 (u r - -u) 
+ 2i0c u sin 3d - u + c 66 (^rr + _ v r - — v J + cn — u e 
+ (co 2 - /3 2 c 44 )i; + 2i(3c u cos 30 - u fl + 2i/3c 14 sin 30 /u r - - u) 
- i^(ci 3 + c 44 ) - wo + C14 cos 30 (w rr --w r -— Woo J 
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- 2ch sin 35 (- w r6 --^w 9 j - i^r(i x5 + i zi ) - $ e ~ re y2 
X cos 35 ($ rr # r $ ee ) + 2re y2 sin 36 (- $ r6 — - * t \ = 0, 

- l"|8(Ci3 + c 44 ) (u r + - u J + 2c 14 cos 36 f- u r g - — Ug\ 

+ en sin 36 (u rr - — u ee - - u r + — u J 

- ifiitu + c 44 ) - u e + c u cos 36 I v„ - "7 v$$ - - v r + — u ) 

r \ r 2 r r l I 

- 2c 14 sin 35 (- y r<? - — u g J + C44 (w rr + -w r + — w gB ) 

+ (w 2 - 2 c 33 )w + re x5 ($ rr + - * r + - # M ) - r/3 2 e z3 * = 0, 

- W(e x5 + e z i) lu r + -u\ - 2e y2 cos 36 (- u rd - — u» j 

- 2 y2 sin 35 / u rr --^uoo u r + — u j 

- i/?(e x5 + ^ 2 i) - 1>0 - e y2 cos 35 / u rr - — u 90 - - v r + — v J 
+ 2e y2 sin 35 (- v r6 - — u# j + e x5 /u; rr + - «; r + — m; 9< , j 

" (3%3W - e xx (*Vr + ^ <*>r + ^ **) + 2 € 22 <*> = 0, (17) 



where 

e 



t = ~, (18) 

cc 



and u, v, and w are the radial, azimuthal, and longitudinal components 
of the displacement vector u. 
In dimensionless form, the boundary conditions (12) to (15) are 

Free surface: 

Cuu r + 612(1* + u e ) — ifUcuW — i0r£zi$ + cos 35[ci 4 (— i(3v + w e ) 

— Te y 2$e] + sin 36[cu(—ipu + w r ) — Te y2 $ r ] = 0, 

Cqg{u + v r -v) + cos 36[c u (-ifiu + w r ) - re y2 ^ r ] 

- sin 36[c u (-i0v + w e ) - Te y2 ^e] = 0, 
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Cui—i&u + w r ) + re X 5$ r + 6 14 cos 39(ug + v r — v) 

+ C14 sin S6(u r — u — ug) = at r = 1. 

Clamped surface: 

u = v = w = 0atr = l. 

Short-circuit: 

$ = at r = 1. 

Operc-circui£: 

— «**$r + &x5(—i&u + a> r ) — c y 2 cos S6(ug + v r — v) 

+ e y2 sin 36(u - u r + v ) = at r = 1. (19) 

III. PERTURBATIONS AND POTENTIALS 

At any given frequency co, we wish to solve the differential equations 
(17) and boundary conditions (19) for the elastic displacement compo- 
nents u, v, and w, and for the electric potential $; these are functions of 
r and 6. We also need to determine the propagation constant (3. Unfor- 
tunately, we have been unable to obtain an exact solution. We shall find 
an approximate solution by combining two techniques which were ap- 
plied successfully in earlier papers. 7,8 First, we observe that eqs. (17) and 
(19) have an exact solution if cu = e y i = 0. 7 In this case, the crystal is 
a member of the hexagonal 6mm class. We make an infinite series per- 
turbation expansion about any modal solution to that problem. This 
results in systems of differential equations and boundary conditions for 
the perturbation contributions to the elastic displacement and electric 
potential. Second, we write these perturbation contributions in terms 
of certain potential functions. The differential equations then decouple. 
With the aid of the boundary conditions, the potential functions can be 
determined; perturbation contributions to the propagation constant can 
also be found. 

The perturbation technique has been used to describe vibrations of 
a sapphire fiber. 8 The equations describing that crystal can be obtained 
from eqs. (17) and (19) by setting <f> and the components of the piezo- 
electric stress matrix e to zero. 

The potential function technique used here is the same as the one used 
in obtaining an exact description of the vibrations of a fiber in the hex- 
agonal 6mm class. 7 

For lithium niobate, we find from the definition (16) and the numerical 
values for the stiffness coefficients 12 that c i4 * 3.6 X 10~ 2 . We will use 
c 14 as a perturbation parameter. This is reasonable since it is small 
compared to one. Instead of treating e y 2 as a separate perturbation pa- 
rameter, we write it as a constant multiple of C14: 
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e y2 = #14- (20) 

For lithium niobate, it turns out that e y2 » 6.8 X 10" 1 and £ « 18. 12 The 
perturbation scheme would work better if e y2 were smaller than this. It 
effectively is, in three out of four differential equations and in all but the 
open-circuit boundary condition, for it is then multiplied by the di- 
mensionless constant t * 1.4 X 10 -1 . In the remaining differential 
equation and boundary condition, however, e y2 is not multiplied by a 
small constant in this fashion. How rapidly the perturbation series ac- 
tually converges will have to be determined numerically. 

We first make a perturbation expansion for the propagation con- 
stant: 

0- t (c 14 ) m t3 m . (21) 

m=0 

When we make a perturbation expansion for the elastic displacements 
and electric potential, it is convenient also to make a Fourier expansion 
in 0. Because of the three-fold symmetry of the crystal about the z axis, 
the Fourier expansion only needs to include multiples of 30, rather than 
6. With Z used to represent u, v, w, or <1>, we assume that 

Z(r,6)= £ (c u ) m £ e iN0 e i3n0 Z m,n( r ), (22) 

m=0 n=— « 

To begin the perturbation scheme, we choose (for m = n = 0) u 0,0 (r) 
e iN0 , . . . , $°'°(r)e iN0 , and /3q to be a modal solution to the unperturbed 
problem, i.e., that for a hexagonal 6mm crystal. N can be any integer. 
It determines which type of modal solution is being considered. N = 
corresponds to an azimuthally symmetric mode, \N\ = 1 to a flexural 
mode, and |iV| > 1 to a higher-order flexural mode. For m = and n j* 
0, set w°- n , . . . , $ '" to zero. The problem then is to determine u m>n (r), 
. . . , $ m '"(r), and fi m for m > 0. We will see that the displacement and 
electric potential contributions vanish when \n\ > m. The functions thus 
need only be determined in the "triangular" region m = 0, 1, 2, 3, . . . and 
|n| < m. 

We next write the perturbation contributions to the elastic displace- 
ments and electric potential in terms of certain potential functions: 

u m,n {r) = SL £ ^.«(r)--^' n (r), 
dr t=\ r 

lr e=\ dr J 

w m ' n (r) = i E H<^?' n (r), 
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$ m > n (r) = i E ntfT n (r), (23) 

with 

s = 3n + N. (24) 

The ne and t\ ( , £ = 1,2, 3, are constants (independent of m and n) which 
must be determined. 

Substitute the potential functions defined in (23) into the perturbation 
expansions. Substitute these, in turn, into the differential equations (17). 
After considerable algebra, we find that the terms multiplied by 
{c l4 ) m e i{3n+N)e yield the following system of differential equations for 

w n , ..., n- n - 

dr e=i 

+ [(co 2 - 2 o c 44 ) + fai/tfis + C44) + Affjrto* + e 2 i)]^ n ) 

- " [c 66 V^4 m '" + (co 2 - j»ii)#rt = Ff- B (r), (25) 

r 

- E KuVj^r 
r/=l 

+ [(u> 2 - ^c 44 ) + /W^is + C44) + PovtT(i x5 + e zl )Wt- n \ 

~ 7" [c 6 6V.^4 m - n + (CO 2 - /8§«44 W' B ] = Wr), (26) 
dr 

E M44 + rj e re xb - ^ (ci3 + c^V 2 ^" 

+ E [(co 2 - 188*33)/*/ - fffcrtTttW-" = ^"to (27) 
E fee**/ + Me x5 + ««l) " «x5M/]V,V7' B 

+ E Wfe»w + PfawtW = F? n (r), (28) 
/=i 



with 



d 2 Id s 2 
dr 2 r dr r 2 



The functions Ff- n {r), j = 1, . . . , 4 are written in the Appendix. They 
are written in terms of functions which have been determined in earlier 
stages of the iterative procedure. When n = 0, they also involve the 
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constant /3 m , which must be found. The functions vanish when m = 
or when m > and \n\ > m. 

In a similar fashion, the perturbation procedure yields a system of 
boundary conditions. 

Free surface: 

? [ 6n dr~* + 6u (fr " S2 ) + iMl3M ' + 0oT§zlVf ]^' n 

- 2c 6G s (£ - l) ft* = Kf" at r = 1, 

a * [£ 2s (I - o »■ - (£ - 1 + s2 ) *■] - **■ 



E [«440*/ " 0o) + ^e x5 ] j- *y- n + PodusW 1 = Kf- n at r = 1; 
/=i ar 



(30) 



Clamped surface: 



Short-circuit: 



e=\dr 

3 d\l/'T' n 

s L ^• n --5 i - = 0atr = l, 
t=\ dr 

E M/^' n = 0atr = l. (31) 



£ i7^? n = 0atr = l, (32) 

Open-circuit: 

t, HxxV* + e x5 (^ - /Jo)] ^j^ + M x5 sW n = Kf- n at r - 1. 

(33) 

The constants KJ***, y = 1, . . . , 4 are written in the Appendix. Like the 
Ff ,n (r), they vanish when m = and are known when m > 0; when rc = 
0, they also involve (3 m . 



IV. SOLUTION OF THE DIFFERENTIAL EQUATIONS 

We now show how to decouple the differential equations (25) to (28) 
and solve them. First, let 
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+ A*f (*u + C44) + PovMe X 5 + «rt)WL (34) 

WW = 6h^T" + (*> 2 " *44W7' B . (35) 

Then by using (25), (26), and procedures similar to those exhibited in 
Ref. 8, we can show that, except in a certain special case to be discussed 
later, 

HT' n (r) = %r s C r x-°[FT ,n (x) + Fp n (x)] dx 

+ %r~ s C r x»[F?- n (x) - F? n (x)] dx, (36) 
Hf' n (r) = %r s C' x~ s [FT- n (x) + F$- n (x)] dx 

- %r" s C x*[F?- n (x) - Ff n (x)] dx, (37) 

Now consider eqs. (27), (28), and (34). They are equivalent to the three 
decoupled equations 

Vfyf" + pty? ,n = f?' n , / - 1. 2, 3, (38) 

provided that 

p? = [(a) 2 - 188644) + /Wci3 + c 44 ) + fart** + t«i)]/tu, (39) 

from (34), 

- p] [M/C44 + vfre x5 - 0o(cu + C44)] 

+ [(*> 2 " /*8s)M " /?o 2 c 23 r»7/] = 0, (40) 

from (27), and, from (28), 

- Pt [ixxVe + Po(e x5 + e zl ) - e x5 ti(] + [-pohzVi + 0oe 2 3M/] = 0. 

(41) 

These imply that the p 2 satisfy the cubic equation 

HxxP) + flfcJKtilP} + 1*44 ~ ^ 2 )(C44P? + /*33 " <»*) 

- P 2 /3 2 (ci3 + C44) 2 ] + r(e x5 p 2 e + 0o 2 e 23 )[(cuP 2 + Pfcu " " 2 ) 
X (2 xS p 2 , + fiUza) ~ 2p 2 2 (ci3 + C44)(e x5 + e zl )] 

+ rp 2 Me X 5 + e 2l )HcuP 2 + /3gc 33 - " 2 ) = 0, (42) 
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and that 

- P 2 (e*5 + e zl )[Te x5 (e x5 + e zl ) + e XI (ci 3 + c 44 )]} 
X |j8o(^i3 + C44)(e 2 3ixx ~ e x5 c zz ) - re x5 (e x5 + e zl )(e x5 p? + $$ gS ) 

~ 6xx(e,5 + «,l)(«44P? + 1*33 - co 2 )}" 1 (43) 

and 

V/ = [cnP 2 + /3 2 c 44 - co 2 - fiePoicu + c 4 4)]/[PoT(e x5 + e z i)]. (44) 

Furthermore, by eqs. (34), (27), and (28), the ff' n (r), t - 1, 2, 3, must 
satisfy 

*u E /?■" = Hf-", (45) 

£ [/vc 44 + Vere x5 - /3 (c 13 + c 44 )]ff n = Ff", (46) 

£ fe*1/ + 0o(ex5 + e zl ) - e x5 n e ]f? n = W- (47) 

These equations can be solved for the ff ,n (r). By (35), eq. (38) also holds 
when £ = 4, provided that 

p 4 2 =(a ) 2 -^c 44 )/c66, (48) 

f m,n = H ?- n /c 66 . (49) 

The next step is to solve the uncoupled differential equations (38). 
The functions ff ,n {r), j = 1, . . . , 4, are either determined completely (n 
5* 0) or else involve /3 m in a known way (n = 0). Using the fact that \pj l,n 
is bounded at r = to evaluate an integration constant, we have as a 
solution to (38) 

tff B (r) = [*Af n + ^ j^ *y.(p/*)/p"tx) dx] J.(p^) 

+ f f xJ s { Pj x)ff- n {x) dx Y s ( Pj r) if pj > 0, 
tff "(r) = \Af n - f 1 xK a (qjx)ff- n (x) dx\ I s (qjr) 

xl s { qj x)ff- n (x) dx K s ( qj r) if pj ^ -q] < 0. (50) 



s: 



For any values of m and n, there are four constants Af ,n , j = 1, . . . , 
4, which remain to be evaluated. When n = 0, (i m must also be found. In 
the next section, we will show how to apply the boundary conditions to 
evaluate these constants. 
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It appears from (36), (37), (50), and (61) that a double integration must 
be performed computationally to obtain H^ n (r) and H^- n (r). Use of (38) 
and integration by parts, however, can reduce this to a single integra- 
tion. 

There is one special case for which the above analysis is not quite 
correct. By using arguments similar to those in Ref. 8, we can show that 
if p 2 ( = for some t, then 

H?' n (r) = C m - n r s + %r< C" x~ s [FT' n (x) + F^x)] dx 

+ fys C r x s [F?- n (x) - Ff n {x)) dx, 
H?- n (r) = C m - n r s + V s C' x-'[FT n (x) + Fp n (x)] dx 

- %r"« C r x s [F?- n (x) - F?- n (x)] dx, (51) 

if n 5* 0. Here C m,n is a constant which remains to be determined. (When 
every pj is nonzero, C m - n is arbitrary in the sense that changing it merely 
changes the constant by which the entire solution is multiplied.) Also, 
in this special case we have when n = 0, 

H r°(r) = CT ,n - f * Fp°(x) dx, 

H?-°(r) = Ct n ~ § l F?-°(x) dx, (52) 

where Cf '" and C^ 1 " must be determined. Now it can be shown that if 
the fiber is vibrating in the lowest-order torsional mode (with u°'° pro- 
portional to r and u°-° = w°<° = <f>°-° = N = 0), then p? = p\ = 0. For this 
case, the solutions of (38) for; = 1 and 4 are 

^f(r) = ["Af"-^- C l x- s+i f]> n (x)dx~\r s 

l r r 
~T I x' +1 ff' n (x)dxr- a ifri9*0, (53) 
2s Jo 

^f°(r) = C x In xf?> °(x) dx + f ' xff-°(x) dx In r if re = 0. (54) 

In eq. (54), an integration constant has been set to zero because it does 
not affect the final solution. It follows from eqs. (23), (43), (44), and (48), 
that when p\ = p\ = and n?^0, the constants A™ ,n and A™-" appear 
in the displacements and electric potential only in the combination Af'" 
- Af- n . Thus for n ^ 0, the constants to be evaluated are A? ,n - A™- n , 
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AT, A^ n , and C m - n . When n = 0, we must find A%-°, Af>", Cf' , C?-°, 
and /3 m . In the next section, we show how to use the boundary conditions 
to determine these constants. 

V. EVALUATION OF THE BOUNDARY CONDITIONS 

For any pair (m, n), there are four boundary conditions, three of which 
are either eqs. (30) or (31), and the fourth of which is either (32) or (33); 
there are also four unknown constants to be found. When n = 0, p m must 
be determined, too. 

When the solutions Pf >n to the differential equations (38) are substi- 
tuted into the appropriate boundary conditions from eqs. (30) to (33), 
a system of equations results which can be written in matrix form as 

J"A m ' n = V m - ra . (55) 

We will not write down here the specific components of these matrices 
and vectors, although it is straightforward to do so. The important things 
to know are the following: The 4X4 matrix J" involves Bessel functions. 
It depends upon /? > but is known once this is determined. The vector 
A m> " consists of the four unknown constants to be determined. The 
vector V m « n contains known constants: Kf ,n , Bessel functions, integrals 
involving fj l,n . When n = 0, it also contains /3 m linearly. 

Incidentally, from a computational viewpoint, it is never necessary 
to differentiate the functions tff '" numerically, either for substitution 
into the boundary conditions or into the functions listed in the Appendix. 
Equation (38) can be used to eliminate all second derivatives of $f' n with 
respect to r. Differentiation with respect to r of the solutions (50) and 
the use of standard relations between Bessel functions and their deriv- 
atives result in analytical expressions for d\l/f ,n /dr. 

The procedure for solving the differential equations and applying the 
boundary conditions is an iterative one. We start with m = 0. We choose 
a modal solution when n = and set all \f/f n to zero when n ^ 0. The tf fl 
satisfy eq. (38) with fj fi = 0. The boundary conditions for this case 
are 

jo A o,o = o. (56) 

From this we obtain for a nontrivial solution the frequency equation 

detJ° = 0, (57) 

which determines /3 as a function of to. This, of course, is the same as the 
dispersion relation for the hexagonal 6mm case about which we are 
perturbing. 

As we iterate on m, we can see from the equations in the Appendix that 
prji,n a K f' n = whenever \n\ > m. It follows that V m >" = in this case 
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and, since det J" ^ 0, that A m -" = 0. Thus all ^f = whenever \n\> 
m. 

If < | n | < m, V m '" is in general nonzero. Since det J n ?* 0, we can 
immediately obtain A mn from eq. (55). 

If n = 0, the analysis is slightly more complicated. It was explained 
in detail in Ref. 8, so we merely give the results here. To obtain /3 m , re- 
place any column of J° by V m " and set the determinant of the resulting 
matrix to zero. The unknown vector A m0 can be written as 

A' n - = C m A ' + D'"' , (58) 

where D m, ° has three unknown components and DJ 1,0 = for some j for 
which Af° ^ 0. Then the equation 

joi)m,o - o (59) 

can be solved for D m0 . Furthermore, C m is arbitrary in the sense that 
varying it varies the constant by which the full solution is multiplied. 
We set C m = 0. 

In this manner, the functions $? *" can be determined iteratively, 
starting with m = 0. For any given value of m, nontrivial results are ob- 
tained only when \n\ <m. The perturbation contributions to the elastic 
displacements and electric potential are then found from eq. (23). The 
full solution is given by eqs. (21) and (22). 

APPENDIX 

Let 

s = Sn + N, 
s+ = 3(n + 1) + N, 
s- = 3(n - 1) + N, (60) 

where n and N are integers. Then 

FT ,n (r) + pm,n {r) m (± _ £^ "j? {^.^ [ £ f$» + *'/»] 

- Pm-j L I(ci3 + c 44 W + r(e x5 + e zl )ve] Vt 

+ L (r^-M^r 1>n+1 L 

e=\ J 
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FT •"(/•) - 



/d s\ m-i r 3 . -| 

Ff n (r) = (- + ") L 7m-,<?44 [ E W - «"J 

3 
- Pm-j E [(Cl3 + C44W + ^(^5 + &zl)Vt] Vt 

*['!-^i-7)]i:5'»-»[i.»*"*«"] 



-iid^-^K- 



2(1 -s+) /d , s+ 



(1^)]} 

[3 "I m-1 T „ 3 

E ^- 1 ' n+1 " ^r M+1 + E ft»-,-(*13 + £ 44 )v s 2 E W 
^=i J y=o L e=\ 

+ 7m-; E (C33M^ + ^z-mtWt - 

[3 "I m-1 f 3 

E W 1Ml - tr 1M1 \+ E A»-,-<«,6 + e 2 i)v s 2 E W 1 
e=i J y=o L ^=i 



+ y m -J E (e 23 ^ - WeWA (61) 
/=i J 



where 



7m = E PjPm-j- 

>=0 



(62) 



^-ig^w((|-.-)[j:^+^] 

-(l + -)LI^ +i -^ +i ]i 
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-^-^>[(s-«-K u - , -(£ + «o*r u "] 



m-\ 3 

- E Pm-j L (ci3M/ + ez\ri]eWi n atr = 1, 



+ d + ^)[i* fI -^ 1 ]) 

+ (— + •+) W u * 1 ]atr-l l 

#3"'" = C 44 "E ft.-,/ { Z \sW + (f r ~ s ) &A ~ *M n \ 

+ 1 [vf. - 2d + s _) (| - S _) ] [ £ *,™ + *r->] 

" 2 [ V ' + " 2(1 " S+) (I + S+ ) ] [A * ? " U+1 " ^- 1 ' B+1 ] at r = 1, 
W" = e x5 "e ftn-> { L [«^ B + (J: - *) W"l - «*4"} 

+ %* [v s 2 _ - 2(1 + s _) (^ -•-)][ z ^ _1,n " 1 + *r u - 1 ] 
- u [ vl - 2d - s + ) (j: +«+)][ £ *r 1,B+1 - +r hn+1 ] 



atr-1. (63) 
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